## 在这个项目中主要存在4个可以直接优化点，分别是：
# 1、迭代次数 N_sim，如果修改 N_sim 则需要相应的改变图的位置。
# 2、fitness, 表达式或者系数都可以调整。
# 3、突变比例
set.seed(1)

getwd()

rm(list = ls())
# install.packages('ggplot2')

library(mgcv) 
library(ggplot2)
library(reshape2) 
library(tidyverse)
source("toolbox.R")


######################################
### settings for simulation number ###
######################################

# parameter
mech <- 5 # choose 1,2,3,4 for other mechanism
N_sim <- 8000  # 迭代次数
a <- 0.3
b <- 0.7
alpha <- 0
beta <- 0


# 描述不同fitness定义方式
# for antagonist
F_antagonist <- function(mech){
  E1 <- (1-H_V_A) # n=1 for A_min & P_max
  E2 <- (1-H_V_A)-alpha*(1-H_V_Po)
  E3 <- (1-H_V_A)
  E4 <- (1-H_V_A)+beta*(1-H_V_Po)
  E5 <- (1-H_V_A)+beta*(1-H_V_Po)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for pollinators
F_pollinators <- function(mech){
  E1 <- (1-H_V_Po)   # n=1 for A_min & P_max
  E2 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E3 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E4 <- (1-H_V_Po)-alpha*(1-H_V_A)
  E5 <- (1-H_V_Po)+alpha*(1-H_V_A)
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

# for plant
F_plant <- function(mech){
  E1 <- (a*H_A_V+b*(1-H_Po_V))
  E2 <- (a*H_A_V+b*(1-H_Po_V))
  E3 <- (a*H_A_V+b*(1-H_Po_V))
  E4 <- (a*H_A_V+b*(1-H_Po_V))
  E5 <- (a*H_A_V+b*(1-H_Po_V))
  
  E <- list(E1,E2,E3,E4,E5)
  return(E[[mech]])
}

###################################################
## load field observation matrix AP_obs, PV_obs ###
###################################################

# # data 应为实际观测值
# AP_obs  <- matrix(data = NA,nrow = 48,ncol = 11) # 这里设置植食昆虫有50种，榕树45种，有机化合物244种
# PoP_obs <- matrix(data = NA,nrow = 18,ncol = 11)
# PV_obs  <- matrix(data = NA,nrow = 11,ncol = 22)

AP_obs  <- read.csv("data/AP_Africa_B.csv", header = TRUE, as.is = TRUE, row.names = 1)
PoP_obs <- read.csv("data/PoP_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)
PV_obs  <- read.csv("data/PV_Africa.csv", header = TRUE, as.is = TRUE, row.names = 1)


############################
### Preparation:Sample #####
############################

index_choose <- c(35,40,16,40,74,12,14,15,1,19) # 按照原则一抽取的 n 种 VOCs
PV_obs_1 <- PV_obs[,c(index_choose)]
index_all <- 1:ncol(PV_obs)

PV_obs_2 <- PV_obs[,c(index_all[-index_choose])]
PV_obs_3 <- PV_obs_2[,c(colSums(PV_obs[,index_all[-index_choose]])<=3)]

# index_1 <- 1:ncol(PV_obs_3)
# index_ <- sample(index_1 , 56)
# print(index_)
# PV_obs_4 <- PV_obs_3[,c(index_ )]

# PV_obs <- cbind(PV_obs_1,PV_obs_4)
PV_obs <- cbind(PV_obs_1,PV_obs_3)


########################################

nA  <- nrow(AP_obs) # 植食动物的数量 （行计数）
nP  <- ncol(AP_obs) # 植物的数量 （列计数）
nV  <- ncol(PV_obs) #  VOC 的数量
nPo <- nrow(PoP_obs) # 传粉蜂的数量

## 生成模拟的初始化矩阵
# 1、初始矩阵不影响模拟结果
# 2、初始时要保证任意植物至少有与一种植食动物有关；任意VOC都至少由一种植物释放。因此下述的colSum 的所以值都大于 0
AP  <- matrix(rbinom(nA*nP,1,0.5),nA,nP)
PV  <- matrix(rbinom(nP*nV,1,0.5),nP,nV)
PoP <- matrix(rbinom(nPo*nP,1,0.5),nPo,nP)
colSums(AP) #  确保 colSum 的所以值都大于 0
colSums(PV) 
colSums(PoP)
# 每一次循环基因变异的比例，该值仅影响收敛速度，不影响结果。
M1 <- 0.2*nP*nV/2 # PV 矩阵突变的比例
M2 <- 0.2*nA*nP # AP 矩阵突变的比例
M3 <- 0.2*nPo*nP # PoP 矩阵突变的比例

#######################################
### Analyses1: Simulation process #####
#######################################
# create plot
par(mar=c(4,4,2,9))
plot(-1, xlim = c(0,N_sim), ylim = c(0,1), ylab = "Fitness", xlab = "Time")
legend(x=N_sim+N_sim/50, y=0.4, title = "Simulated", legend = c("Plant Fitness","Antagonist Fitness","pollinators Fitness"),
       pch = c(3,4,5), col = c("red","red","red"), xpd=TRUE, bty="n",title.font = 2) 
legend(x=N_sim+N_sim/50, y=1, title = "Parameter", legend = c(paste0("a = ",a), paste0("b = ",b), paste0("alpha = ",alpha), paste0("beta = ",beta)),
       xpd=TRUE, bty="n",title.font = 2) 


# 根据实地观察，计算并绘制熵值和fitness
{
  AV_obs  <- as.matrix(AP_obs) %*% as.matrix(PV_obs)
  PoV_obs <- as.matrix(PoP_obs) %*% as.matrix(PV_obs)
  # "H_A_VOC" 函数 i:矩阵 o:互信息/条件熵，toolbox
  # H_V   <- H_A_VOC(PV_obs)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV_obs)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV_obs)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV_obs)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV_obs)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP_obs)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP_obs)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV_obs)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV_obs)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP_obs)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP_obs)[["Hn_V_S"]]

  # 对抗系统中的 fitness 表达, toolbox
  E_plant       <- F_plant(mech)
  E_antagonist  <- F_antagonist(mech)
  E_pollinators <- F_pollinators(mech)



  # E_obs <- c(0, E_plant, E_antagonist, E_pollinators, H_V,
  #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  # names(E_obs) <- c("N","E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
  #
  ## 野外实际观测值划线，图A
  # 条件熵 for PV, AV and AP
  # abline(h = H_P_V, lty = "dotted", pch = 8)
  # abline(h = H_A_V, lty = "solid", pch = 8)
  # abline(h = H_A_P, lty = "dashed", pch = 8)
  ## 植物和植食动物的fitness，图B
  abline(h = E_antagonist, col = "red")
  abline(h = E_pollinators, col = "blue")
  abline(h = E_plant, col = "green3") # the same as H(A|V) in mech = 1
}


# 模拟过程
A_PV  <- list()  # 每次模拟后存储 PV矩阵 的 list
A_AP  <- list()  # 每次模拟后存储 AP矩阵 的 list
A_PoP <- list()  # 每次模拟后存储 PoP矩阵 的 list
# E = matrix(NA, nrow = N_sim +1, ncol = 15) # E 用来存储每个模拟之后的变量值 N*14
# colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators", "H_V", "H_P_V", "H_V_P", "H_A_V", "H_V_A", "H_A_P", "H_P_A","H_Po_P","H_P_Po","H_Po_V","H_V_Po")
E = matrix(NA, nrow = N_sim +1, ncol = 8) # E 用来存储每个模拟之后的变量值 N*14
colnames(E) <- c("N", "E_plant", "E_antagonist", "E_pollinators",  "H_A_V", "H_V_A","H_Po_V","H_V_Po")


## 基于随机初始化的第一次模拟值
n = 1
{
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # "H_A_VOC" 在toolbax.R 中输入矩阵，输出互信息和条件熵
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # 绘图
  # points(n,  H_P_V, col = "red", pch = 24) 
  # points(n,  H_A_V, col = "red", pch=  21)
  # points(n,  H_A_P, col = "red", pch=  22) 
  # 图B
  points(n,  E_pollinators, col = "blue", pch = 5)
  points(n,  E_antagonist, col = "red", pch = 4)
  points(n,  E_plant, col = "green3", pch = 3)
  
  #写入 E, A_AP, A_PV 的第N次
  # E[n,] <- c(n, E_plant, E_antagonist, E_pollinators, H_V,
  #            H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[n,] <- c(n, E_plant, E_antagonist, E_pollinators,
             H_A_V, H_V_A, H_Po_V, H_V_Po)
  
  A_AP[[n]] <- AP
  A_PV[[n]] <- PV
  A_PoP[[n]] <- PoP
}


### 植食动物和植物依次使其 fitness 最大化 迭代模拟
for (n in 1:(N_sim/4)) # N_sim : 模拟次数。由于模拟过程包括植食动物、传粉蜂的基因突变和植物两次基因突变四个部分，因此需要 "/4"
{
  print(n)
  ## 传粉蜂：最大化fitness
  for (m in 1:M3) {
    PoP_new <- random.sample_1element1(PoP) #function details see toolbox.r
    PoV <- PoP_new %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP_new)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP_new)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_pollinators_new > E_pollinators)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PoP <- PoP_new 
    }
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  
  E[4*n-2,] <- c(4*n-2, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  
  A_AP[[4*n-2]] <- AP
  A_PV[[4*n-2]] <- PV
  A_PoP[[4*n-2]] <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  ## 被子植物：最大化fitness
  for (m in 1:M1) {
    PV_new <- random.sample_1element1(PV) #function details see toolbox.r
    PoV <- PoP %*% PV_new
    AV  <- AP %*% PV_new
    # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_plant_new > E_plant)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PV <- PV_new
    }
    
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n-1,] <- c(4*n-1, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[4*n-1]] <- AP
  A_PV[[4*n-1]] <- PV
  A_PoP[[4*n-1]] <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  
  ## 植食动物：最大化fitness
  for (m in 1:M2) {
    AP_new <- random.sample_1element1(AP) #function details see toolbox.r
    AV <- AP_new %*% PV
    # H_V   <- H_A_VOC(PV)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP_new)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP_new)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
    
    # fitness function
    E_plant_new       <- F_plant(mech) 
    E_antagonist_new  <- F_antagonist(mech) 
    E_pollinators_new <- F_pollinators(mech) 
    
    # 基于优化机制判断是否保留此次基因突变
    if(E_antagonist_new > E_antagonist)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    AP <- AP_new 
    }
  }
  
  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]
  
  # fitness function
  E_plant       <- F_plant(mech) 
  E_antagonist  <- F_antagonist(mech) 
  E_pollinators <- F_pollinators(mech) 
  
  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[(4*n),] <- c((4*n), E_plant, E_antagonist, E_pollinators, H_V, 
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n,] <- c(4*n, E_plant, E_antagonist, E_pollinators,
               H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[(4*n)]] <- AP
  A_PV[[(4*n)]] <- PV
  A_PoP[[4*n]]  <- PoP
  
  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24) 
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonist, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
  
  ## 被子植物：最大化fitness
  for (m in 1:M1) {
    PV_new <- random.sample_1element1(PV) #function details see toolbox.r
    PoV <- PoP %*% PV_new
    AV  <- AP %*% PV_new
    # H_V   <- H_A_VOC(PV_new)[["Hn_V"]]
    # H_P_V <- H_A_VOC(PV_new)[["Hn_S_V"]]
    # H_V_P <- H_A_VOC(PV_new)[["Hn_V_S"]]
    H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
    H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
    # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
    # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
    H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
    H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
    # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
    # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]


    # fitness function
    E_plant_new       <- F_plant(mech)
    E_antagonist_new  <- F_antagonist(mech)
    E_pollinators_new <- F_pollinators(mech)

    # 基于优化机制判断是否保留此次基因突变
    if(E_plant_new > E_plant)
    {E_antagonis  <- E_antagonist_new
    E_plant       <- E_plant_new
    E_pollinators <- E_pollinators_new
    PV <- PV_new
    }

  }

  # 更新矩阵，并计算互信息量和条件熵
  AV  <- AP %*% PV
  PoV <- PoP %*% PV
  # H_V   <- H_A_VOC(PV)[["Hn_V"]]
  # H_P_V <- H_A_VOC(PV)[["Hn_S_V"]]
  # H_V_P <- H_A_VOC(PV)[["Hn_V_S"]]
  H_A_V <- H_A_VOC(AV)[["Hn_S_V"]]
  H_V_A <- H_A_VOC(AV)[["Hn_V_S"]]
  # H_A_P <- H_A_VOC(AP)[["Hn_S_V"]]
  # H_P_A <- H_A_VOC(AP)[["Hn_V_S"]]
  H_Po_V <- H_A_VOC(PoV)[["Hn_S_V"]]
  H_V_Po <- H_A_VOC(PoV)[["Hn_V_S"]]
  # H_Po_P <- H_A_VOC(PoP)[["Hn_S_V"]]
  # H_P_Po <- H_A_VOC(PoP)[["Hn_V_S"]]

  # fitness function
  E_plant       <- F_plant(mech)
  E_antagonist  <- F_antagonist(mech)
  E_pollinators <- F_pollinators(mech)

  # # 写入模拟完植食动物的 E, A_AP, A_PV
  # E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators, H_V,
  #                H_P_V, H_V_P, H_A_V, H_V_A, H_A_P, H_P_A,H_Po_P,H_P_Po,H_Po_V,H_V_Po)
  E[4*n+1,] <- c(4*n+1, E_plant, E_antagonist, E_pollinators,
                 H_A_V, H_V_A, H_Po_V, H_V_Po)
  A_AP[[4*n+1]] <- AP
  A_PV[[4*n+1]] <- PV
  A_PoP[[4*n+1]] <- PoP

  # # 绘图
  # points((4*n),  H_P_V, col = "red", pch = 24)
  # points((4*n),  H_A_V, col = "red", pch=  21)
  # points((4*n),  H_A_P, col = "red", pch=  22)
  # 图B
  points((4*n),  E_pollinators, col = "blue", pch = 5)
  points((4*n),  E_antagonis, col = "red", pch = 4)
  points((4*n),  E_plant, col = "green3", pch = 3)
  
}





####################
# 互惠对抗1.4版
# 优化代码计算量